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The load carried by a queuing system under equilibrium conditions is 
the average amount of server usage per unit of time. In telephony, this 
parameter is often evaluated by recording the number of busy servers at 
regular time intervals; these readings are then cumulated and their sum, 
after division by the number of observations, is an unbiased estimate of 
the carried load. The purpose of this paper is to derive exact formulas 
for the computation of the variance of this measurement in systems with 
arbitrary input and departure rates. The results obtained here thus apply 
to a wide class of teletraffic models which includes, in particular, the delay- 
and-loss systems with finite- or infinite-source inputs, exponential service 
times, and arbitrary defection rates from the queue. Problems related to 
computations are also considered, special attention being paid to the 
reduction of both computer time and storage when the number of states is 
large. 

I. INTRODUCTION 

Analysis of the stochastic behavior of traffic measurements is of 
considerable practical relevance, as it provides means for appraising 
field data as well as guidelines for selecting performance standards. 
Load measurements play a central role in this effort, and determina- 
tion of their accuracy is therefore of particular interest. The present 
investigation yields an answer to this problem for a broad class of 
teletraffic models. 

Whenever statistical equilibrium prevails (and it is assumed to 
throughout this paper), the load carried by a service system is the 
average amount of server usage per unit of time or, equivalently, the 
average number of busy servers at an arbitrary instant. In telephony, 
an estimate of this parameter is often obtained by "switch-counting." ■ 
This statistic, which is determined by recording the number of busy 
servers at regular intervals and then by taking the arithmetic mean 
of these discrete observations, is an unbiased estimate of the carried 
load. 
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The variance of this measurement, called hereafter the switch-count 
load to distinguish it from the estimate obtained by continuous observa- 
tion, was first determined approximately by Palm 2 and Hayward 1 in 
the case of an infinite server group with Poisson input and exponential 
holding times. This result was later extended by Benes, 3 who obtained 
the exact variance of the switch-count load for groups of finite sizes 
without waiting positions (loss systems). A further generalization to 
loss systems with recurrent input and exponential service is due to 
Neal and Kuczura. 4 Their formal analysis stops, however, with a 
derivation of the Laplace transform of the covariance function of the 
underlying carried-load process. From this point on, they proceed 
numerically, since explicit inversion of the transform appears to be 
difficult in general. 

In this paper we are concerned with derivations of exact formulas 
for the variance of the switch-count load in finite systems with arbitrary 
state-dependent input and departure rates. The results presented here, 
therefore, fill a rather large gap, since they apply to a broad class of 
teletrafhc models that includes, in particular, the (finite) delay systems 
with exponential holding-time distributions, arbitrary defection rates 
from the queue (if one is allowed to form) and either Poisson or 
quasi-random input (in the latter case, the traffic is generated by a 
finite number of sources that place demands for service at the same 
constant probability rate when free but that do not submit requests 
while being either served or waiting). 

Let N(t), the state of the system at time t, be defined as the number 
of busy devices at that instant (by device, we mean here either a server 
or a waiting position). Let c and d be, respectively, the number of 
servers and the number of devices. 

Unless stated otherwise, we make the following assumptions : 

(i) When N(t) = n and ^ n < d, the probability that a re- 
quest originates during (/, t + h), h > 0, is of the form \„h 
+o(h), with X„ > 0. 

(ii) The requests which are submitted when all the devices are 
occupied are dismissed and, accordingly, \ d is set equal to zero. 

(in) When N(t) = n and < n ^ c, the probability that a service 
time terminates during (t, t + h) is of the form n„h + o(h), 
where n„ > 0. 

(iv) When N(t) = n > c, the probability that either a service time 
terminates or a waiting request defects from the queue is of 
the form n n h + o(h) where p„ > and n ^ d. 
(v) When a server becomes free, it is immediately reseized by one 
of the waiting requests if any are present in the system at that 
time. 
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Let N c (t) be the number of busy servers at time t and let c a n be 
the smaller of the two integers a and n. Then 



N c (t) = c a N(t) = ] 



N(t) if N(t) ^ c, 
c if N(t) > c, 



and the switch-count load, L n (T), based on n observations (scans) 
made over [0, T~\ at times t, 2r, • • • , nr, is, by definition, equal to 

"- 1 t N c (jr), 

where t = T/n. 

Let Cov [A^ c (ii), N c {to)~\ be the covariance between N c (ti) and 
N c {t2). Under equilibrium conditions, this covariance depends only on 
1 1\ — h | so that 

Cov [iV«(<x), N c (t 2 n = Cov [iV c (0), iST c (|<i - «.!)]. 

Hence, the variance of L„{t), cast in a form that will be convenient 
later, is given by the formula (Ref. 3, p. 137) : 



VarL„(T) = n" 2 E (n - \k\)R e (kr), (1) 



fc— n 



where 



fl £ (A;r) ■ Cov [iVe(O), A^ c (A;t)] 

= Cov[A^ c (0), AT.(|*|t)]. 



It is clear from (1) that the variance of the switch-count load is 
completely determined by the covariance function R c (-) of the carried- 
load process [N c (t), —<*><t<oo] f and therefore much of what 
follows is concerned with expressing R c (-) in the most convenient form. 

The covariance function can be stated at first in terms of the transi- 
tion probabilities, and the resulting expression can then be reduced 
by taking the structural properties of the process into account. But 
alternate forms can also be obtained by making use of the fact that 
the conditional expectations, E{N e (t)\N(0) = m], m = 0, 1, •••, d, 
satisfy simple linear differential equations. The covariance formulas 
obtained by these diverse procedures exhibit distinct features that may 
be exploited in the computations. In all cases, however, R e (t) is ex- 
pressed as a diagonal, positive-definite quadratic form which reveals 
that R c (-) is completely monotonic. 5 

Expressions for the transition probabilities, the covariance function, 
and the variance of the switch-count load are derived in Sections II, 
III, and IV, respectively. The variance of load measurements based on 
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continuous observations can be found in Section IV. Extensions of 
the results of Sections II to IV to reversible Markov processes are also 
considered. Questions of a computational nature are dealt with in 
Section V, while Section VI is concerned with some asymptotic 
properties (d large) of the spectrum of the underlying process. 

The formulas presented here have been programmed and used to 
explore the effects of parameter changes on the variance of the switch- 
count load. The result of that investigation will appear in another 
paper. 

II. TRANSITION PROBABILITIES 

In this section, we express the transition-probability matrix as a 
symmetric product of vectors and matrices. As becomes apparent later, 
this representation makes it possible to write the variance of the switch- 
count load in a way that greatly simplifies its evaluation. 

Let p mn (t) be the probability of a transition from state m to state n 
in time t : 

Vmn (t) m Pr [N(t) = n\N(0) = to], m, n - 0, 1, ■ ■ ■ , d. 

These transition probabilities satisfy the following system of 
differential equations : 

d 

-T. PmO = MlPml — AoPmO, 

jfPi** = Mn+lpm,n+l + K-lPm,n-l — (X« + Pn)Pmn, 1 = n < &, 

-jj Pmd = Xd-lPm.d-1 — fJ-dPrnd- (2) 

Let 

— Ao Ao 

Mi — (X1+M1) Ai 

pi — (X2+/i2) X2 



A rf = 



— (Xd-1+Md-l) Xd_l 

Hd —pa J 



and 



p d (t) s [p mn (0], m, n = 0, 1, -,d, 



be the transition-probability matrix. Capital and lower-case bold-face 
letters are used exclusively to designate matrices and vectors, re- 
spectively. 
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With this notation, the system of differential equations (2) becomes 

j t P d (t) =P d (t)-A„, t^O, (3) 

so that, for A: = 1,2, • • •, 

^ k P d (t) = ^l i P d (t)-A d , tZO. (4) 

It follows from our assumptions that if the system is in state m at 
time zero [A r (0) = w], then \\m t ioPmm(t) = 1 and limao/> mn (0 = 
for n 7^ m. Hence, with Id the identity matrix of order d + 1, the 
initial conditions take the following form : 

P rf (0) ^ limP rf (0 = I,. 
and by (3) and (4) we therefore have 



lim^P d (0 = AS. (5) 

no 



dt k 



The initial conditions state that Pd(-) is right-continuous at t = and 
imply that P<*( •) is continuous for all t > 0. By (3) and (4), all the de- 
rivatives of Pd(-) exist for t > 0, and by (5) they are also right- 
continuous at t = 0. An application of Taylor's theorem then yields 
(Ref. 6, pp. 240 ff.) 

P rf (0 = exp (At) = t h A &> ' = °- ( 6 ) 

The elements of A d situated immediately either above or below the 
diagonal are all strictly positive and so Ad can be symmetrized. In- 
deed, let 

D rf = diag [5 , Si, • • • , 5 d ] 

with 

5o = r and 8 m = f ( /' Ma " \ M "' ) = S?>m*, m = l,---,d, 

\ AoAi ■ • • A m _i / 

where (i) f is a nonvanishing but otherwise arbitrary constant, (ii) 
the p m are the equilibrium state probabilities, and (in) £ = fp&. 
Without loss of generality, we can — and shall — set f = pi so that 
£ = 1 and 

D, 1 = diagOfc, p\, ■-,pl]. (7) 

It is easy to verify that 

S,y = D^A^Dd (8) 
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is symmetric, its nonvanishing elements being 

s mm = - (Xm + Mm), m = 0, 1, ■ ■ ■, d, (Ad = 0), 
s m , m+ i = s m+Um = (\ m fim+i) 1 , m = 0, 1, • • •, d - 1. 

Hence, by (8) we have 

A5 = Drf-SJ-D,- 1 , k = 0, 1,2, ••■, (9) 

and, by (6), 

P d «) = Drf-exp (S^-D,- 1 = E h (Dd-SS-D," 1 )^. (10) 

fc=0 # ! 

The representation of k d in terms of the symmetric tridiagonal matrix 
Sd entails substantial formal simplification of the final results. And it 
is also particularly convenient computationally, since the determina- 
tion of the characteristic values of A d (which are needed for an exact 
solution) is best carried out after symmetrization. 

The matrices Ad and Sd clearly have the same characteristic values, 
ro, J'i, • • • , r d . But Sd is symmetric and is therefore unitarily similar to 
the diagonal matrix 

C d = diag [ro, /'i, ■ • • , fi}. 

This means that an orthogonal matrix B<f exists such that 

S d = B'a-Cd-Bd, B d -B d = B' d -Bd = I d , (11) 

where B' d is the transpose of B d . 

But Sd is also tridiagonal, and its off-diagonal elements never vanish. 
Hence, S d is nonderogatory and its characteristic values are necessarily 
distinct (Ref. 7, p. 26). The elements in the ?ith column of B' d are then 
the components of the (uniquely defined) normalized characteristic 
vector associated with the nth. characteristic value r„(n = 0, 1, • • •, d). 

We now substitute (11) into (10). This yields 

P*(0 - t h (Dd-B'd-CS-Bd-Dd- 1 )^ 

so that 

Pd(t) = Dd-B rf -exp (Cd-0-Bd-D d -' 

= Dd-B'd-diag [V<", e'", • • •, e**].B 4 -D?\ (12) 

We note now that all the row sums of Ad vanish and one of the 
characteristic roots, r , say, must therefore be equal to zero. Further- 
more, known extremal properties of the characteristic values can be 
used to show that r\ t r 2 , • • • , r d are negative. It is also readily seen that 

Pd*>= ivlvl •••,Pd) / 

is the characteristic vector of S d that corresponds to the vanishing 
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characteristic root r . Indeed, let e<* and Od be the (d + 1) dimensional 
(column) vectors whose components are all equal to 1 and 0, respec- 
tively. Then, since A rf -ed = Od, we have, by (8) and (7), 

D d Sd D^ 1 e</ = Dd-Sd-Pd*' = Od. 
But none of the diagonal elements of Da vanishes and the relation 

D d -Sd-Pd J) = Od 

can hold if and only if Sd-pd^ = o d . Thus, pp is the characteristic 
vector associated with r ( = 0), a fact that may be of relevance in the 
computations, as a comparison of p£ J) with Dj" 1 provides an accuracy 
check for the method used to determine the characteristic vectors. 

In the derivation of formula (12), advantage was taken of the fact 
that the transition-rate matrix Ad is symmetrizable. It is worth 
noting that this relatively simple expression for Pd is a consequence of 
this property, and therefore holds for all (and actually only for) 
reversible Markovian processes with finite state spaces. Indeed, by 
definition, the class of these processes — which includes those of the 
birth-and-death type — is fully characterized by the following condi- 
tions (R.efs. 8 and 9) : 

PmPm„(0 = P„Pnm(t), VI, 71 = 0, 1, ■ ■ • , d, (13) 

or, equivalently, by the single relation : 

D7 a -P< = P.-D7 2 . (14) 

Hence (12), written in terms of Sd, implies that 

D d " 2 -Pd = Dj-'-exp (S d t) -D* 1 

= (D^-exp (SdO-D," 1 )' = Pd-D*- 2 , 

and (14) is therefore satisfied. 

Conversely, we show next that (14) is a sufficient condition for 
(12) to hold. 

Pre- and post-multiplication of (14) by D d yield 

D7 & -P«-Dj- Dd-Pd-Dd" 1 . (15) 

Substituting the expansion of Pd as given by (6) into (15), and perform- 
ing the multiplications by Dj and D7 1 under the summation sign 
(which is clearly legitimate), we obtain: 

t n (Dd-'-Aj-Dd)^ = f n (Dj-CAd'-Dr 1 )**. t ^ o. 

However, this relation cannot be satisfied unless 
D.r'Ad D,, = Dd-A'd-Dd -1 , 

LOAD MEASUREMENT VARIANCE 1283 



so that, by transposition, 

(Dr'-Arf-D*)' = Dj-'-Atf-D* 

This means that A d is symmetrizable by pre- and post-multiplication 
by Dr 1 and D d , respectively, and (12) then follows as shown earlier. 

Under the assumption that the process is reversible and that all the 
states communicate with each other 10 (i.e., p mn (t) > 0, m, n = 0, • • •, 
d, t > 0), the characteristic roots of A d are necessarily simple. (Note 
that A d , and hence S d = D^-A^D*, need no longer be tridiagonal.) 
This can be proved as follows. 

The matrix S d is symmetric and can therefore be tridiagonalized by 
a method from Householder (Ref. 7, pp. 152, 153, 290-293, and 343). Ac- 
cording to this procedure, the tridiagonalization of S d is achieved by 
successive right and left multiplications by symmetric orthogonal 
matrices, Ui, U 2 , • • • , Ud-i, of the form 

V T = Id - 2w r w;, 

where w r is a suitably chosen d + 1 dimensional (column) vector 
whose first r components are zero. (All the U r are of order (d + 1) 
and U? = Id, r = 1, 2, ■ • •, d — 1.) A direct application of the results 
derived in Ref. 7, above, shows that S rf admits of the following repre- 
sentation : 

Sd = U x -U 2 Ud-i-Td-Ud_i U 2 -Ui, 

where Td is a symmetric tridiagonal matrix of order d + 1 . 

Let On be the elements of Td. 

We are now faced with two possibilities. Either 0,. 1+ i = 0<+i,« ^ 
for i = 0, 1, • • •, d — 1, or there is an index j(<d) such that Bjj+i 
= dj + u = so that 

T y 



Td = 



Td-y-x , 



(1(5) 



In the first instance, all the characteristic roots of S d , and hence of 
Ad, are distinct (Ref. 7, p. 26). To complete the proof, it is therefore 
sufficient to show that the second contingency cannot occur when all 
states communicate with each other. To this end, we proceed in- 
directly. We assume that (16) is satisfied for some j < d and show 
that some states then do not communicate with others. 

When (16) holds for j < d, we have, for any k ^ 0, 



A5 = Dd-Sd-Dd" 1 



= Dd-Ur 



■Ud_i 



m 



Td-,-i 



•Ud-i- 



•Ui-Dr 1 . (17) 
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The first row in D d is (po*, 0, • • •, 0) and the elements in the first 
column and first row of each of the U r 's are zero except for their first 
component, which is always equal to 1. Hence, 

(1,0, ■•■ ! 0)-D d -U 1 Ud_i = (po-'.O, •••,0). (18) 

Similarly, since the U r 's are symmetric, we have 

Ud-i UrDr-d.O, ■••,0)' = (p§, 0, ••-,0)'. (19) 

Hence, by (17) to (19), 

Poo(0 = d,o, ■■-.o>( ft 5 l ifi A *''*) (1 ' ' '"' 0) ' 

= t n ^* W 

fc-o* ! 

where 0o*° is the element belonging to the first row and first column 
of T$. 

Let a mn and s m „, m, n = 0, 1, • • •, d, be the elements of Ad and Sd, 
respectively. Under the present assumptions, a m „-a„ m gj: and 
s m „ = (a m „-a„ m )*, m, ra = 0, 1, ■•-, d, m ^ n. The elements in the 
first r rows and columns of Sd are therefore uniquely determined by the 
elements in the first r rows and columns of Ad. Similarly, the vector w r 
depends only on the components, s mn , of Sd for which either m £| r — 1 
and it = r — I, ■ ■ •, d or, by symmetry, n ^ r — 1 and w = r — 1, 
• ■ •, d (Ref. 7, pp. 290 ff). Consequently, the elements of T* (which are 
all obtained after j — 1 steps) depend only on the elements of the 
first j + 1 rows and columns of Ad. This implies that the transition 
probability Poo(t), as given by (20), is independent of the rates a m „, m, 
it > j. However, the process being Markovian, this can only be true 
if Po m (t) = for in > j which means (since, by assumption, j < d) 
that state does not communicate with states j -\- 1, • ■ ■ , d, as was to 
be proved. 

III. COVARIANCE FUNCTION 
3.1 First version 

The covariance function of the carried-load process is, by definition, 

d 

Re(t) = L (c A n)-(c a m)p„p nm (t) - M 2 cX 

m , n=0 
d 

= E (c a n)-{c a m)p„\_p nm (t) - /),„], 

m , n=0 



where 



M el = EX c {0) = £ (c a ,i)p„. 

71=0 
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However, if 7 is an arbitrary constant, the covariance of the process 
[N e (t) + 7, t ^ 0} is also R c . Hence, with the notation 

p„ = (c a n) +7, n = 0, 1, ■•-, d, 

we also have 



Re(t) = £ pnPmPn[pnm(t) — Pm~\- 
m,n=0 



(21) 



Let 



Po Pi 
Po Pi 

Po Pi 



Pd 



and 



G d (0 = \jp nm (t) - pj. 

The matrix P d can be obtained by letting J— ><» in (13). Hence, 

P d - Drf-Bi-diag [1, 0, 0, • • -, OD-Brf-D," 1 

and 

G„(0 = P d (t) - P d 

= D d -B;-diag [0, «nt f . . ., e^Bd-D*" 1 . 

We now introduce two auxiliary row vectors : 

Td = (po, pi, ■ • •, Pd), s' a = (popo, Pipi, • • *i Pdpd)- 

Then the coefficient of e""' in the linear form 

s' d -G d (t)-T d = s d -D<rB d -diag[0, e r1 ', •••, e^ t ~\-B d -'Dt l u 

is the same as the coefficient of e T < 1 in (21), and we may conclude that 

Rc (0 = Sd-Dd-B'd-diag [0, e r »', • • •, e'^-Bd-D^-r,,. 

With the notation 



we have 
so that 



q'd = (popo, •••, PdPd), 

q„ = s'i-Dd and q d = Di" 1 -^, 



(22) 



B.(0 = q'd-B d -diag [0, e'»', ■ • •, e^-B d -q d 
or, alternatively, 

B.(0 = £ ftfr" 1 (23) 

i-l 

with 6, the ith component of the row vector q' d B d . This last expression 
shows that the coefficient of e rit in either (22) or (23) is necessarily 
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nonnegative. Furthermore, since all the r,-'s (?' = 1, • • •, d) are nega- 
tive, we also have 

(-l)"^- k R c (t)^0, tfcO, *-0, 1, .-., 

and R c is therefore completely monotonic over [0, °o). 5 

If we now set y = — c, the last d — c components of q,? are equal 
to zero, and so the determination of R c by means of (22) necessitates 
only the computation of the first c components of the characteristic 
vectors. Formula (22) is therefore often well-suited to the case of delay 
systems. But unless the number of waiting positions exceeds the 
number of servers, greater reduction of the computations can be 
achieved by means of the formulas derived below. 

In the preceding derivation, the p's are independent of the arrival 
and departure rates, and the formulas of this subsection therefore hold 
for arbitrary, reversible, Markov processes with finite state spaces. 
In contrast, the results of the next subsection are restricted to birth- 
and-death processes. 

3.2 Alternative forms 

Multiplying the nth equation in (2) by (c a n) and then summing 
with respect to w(0 ^ n ^ d), we obtain, after rearranging and 
canceling terms, 



But 



E (C A n)-j Vmn = Z^nPmn ~ t Hn Vm ,, (24) 

;i=0 «t 1=0 N = l 



E (c a n)p mn (t) = E{N c (t)\N(0) = m], 



so that, by (24), 

~E{N e (t)\N(0) = m] = 8 EXnP«- " t PnVmn. (25) 

at «=0 n = l 

Adding and subtracting kE{Nc(J) \N(0) = »/} on the right-hand 
side of (25), we obtain 

jE{N c (t)\N(0) = m) = K E{N c (t)\N(0) = m) + £ rftop™, 
at n=o 

m-0, 1, ■■•,d, (26) 

where 

f (X„ — n„ — K7i) if 11*0,1, • ■ -, c — 1, (mo = 0), 

P«(«) = \ — (kc + He) if » = c, 

[ -kc if n = c + 1, •■•, d. (27) 
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In the preceding formulas, k is an arbitrary real number that may 
be positive, negative, or null. We see later that the covariance formulas 
can occasionally be simplified by appropriate choices of k. 

Taking the initial conditions, 

E{N c (0)\N(0) = m\ = (c a m), ^ m < d, 
into account, the solution of (26) is 
E{N c (t)\N(0) = m) 

= (c a m)e" + /*%««-"> T £ p*( K )p mn (w)|.dM (28) 
so that 

Rc(t) =ZM m)p m E[N e (t)\N(fl) = m\ - M 2 cl 

m=0 

= M c% e" - M 2 cl 

+ f'e«(*-»J I(ca «,)p» T £ P»(*)Pm»(w)l-d«, 

7o m=0 L n=0 J 

where M c2 = EN 2 c (0). By means of (13), the preceding relation can 
be expressed in a much more convenient form : 

R c (t) m M e2 e<t - M 2 cl 

+ /"V«-«> £ P :( K )p„ [f(cA m)p nra (w)l-dM 
= itf c2 e«< - M 2 C1 

+ f e«c«-«> £ p*n(K)p n E{N e (u)\N(0) = n}-d«. (29) 

Next, substituting (28) into (29), we obtain 
R c (t) = M c2 e Kt - M 2 cl + te** E (c a n)-p*(K)-p„ 

+ /" e«c«-»> £ Pm(K)-p„ [ W e'<"-'>pnm(v)>dvdu. 

J0 n,m = ^0 

Let 12* be the Laplace transform of # c and p™„ that of p„ m . The pre- 
ceding relation then yields 

R e (s) = + Tk \2 L ( c A n)-pM-Pn 

S — K S [S — K) z „ = 

+ £ pTW-iCW-P-^^i- (30) 

n,m=0 V. 6 K J 

We know, however, that p nm {t) is of the form 

d 

P»»(0 = Pm + E y nmi e rit , 
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so that 

a d 

Rc(t) = L (c A n)-(c A mj)p» L7»m.e r ''. 

m , n=0 i"=l 

This implies that the only poles of R* c are r„ i == 1, • • •, d, and that 
lim^oo R c (t) = 0. Taking these two facts into consideration, we see 
at once that (30) reduces to 

R*ds) = £ p:(k)- p *M- Pii E 7 '"" ] 



n,m=0 i = l (»i — «) S — J'j 

provided k ^ r,-, i = 1, 2, ■ ■ ■ , d. And referring back to the derivation 
of (22), it is readily seen that 



Rc(t) = [qSOOT-B'a-diag 



[• ■ e Tdt 1 * 



QTlt 

(r~«)" ' ' (r d 

k ^ n, i = 1, •••, d, (31) 

where 

MOOT = OoW-po, ■ •*, pSW-pD- 

The modifications needed when k is equal to one of the characteristic 
roots are immediate. Let 

diag 0) Oo, a h ■■-, a d ~] 

be the diagonal matrix obtained by setting the jth diagonal element of 
diag [a , a h ■ ■ ■ , a rf ] equal to zero. Then if k = r, we must have, with 
some as-yet-undetermined constant a and a% the variance of N c (0), 



, i = 1, •• •, d -Bd-qSCfc). 



B e (0 = (<r 2 c + a)e" 




+[q!W]'-Bi-diag«> 




But R e (0) = <r% so that 





,dj 



a = - [qJWJ-Bi-diagW [0, (r< - k)~ 2 , i = 1, • • •, d]-B d -q2(/c). 
Hence, 
B.(0 = <^«< + [qSWJ-Bi-diagW [o, ^ I ^ , * = 1, ■ ■ ■, d\ 

■Bg.qJW- (32) 

It should be noted that (32) is valid even if k ^ ?-,-, i = 1, • • • , d, 
and that it should be used in the computations [rather than (31)]] 
whenever k is "close" to one of the characteristic roots so as to avoid 
accuracy losses (see below). Since k is arbitrary, one could always 
choose it so that it is not "close" to any of the characteristic roots. 
But, as shown next, it is often preferable to select it in such a way as 
to reduce the amount of computation, and this, in turn, may dictate 
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(33) 



the use of (32). As we have seen, 

Rc(t)= [qSOOI-B^diagfo, fo'Z'^ i ••-! (rZ-c)' ]' 8 *'^^- 

But, as noted previously, p^ }) is the characteristic vector associated 
with the vanishing root (?•<>), so that 

(W-K-dba [o, pj^j, , .... ^^J-B.-pi*. = o. 

Consequently, (31) remains valid for all [qSOO]' of the form 

{[poOO + t]-Po, • ■ -i La* 00 + 7]-7>1}i 

where y is an arbitrary constant. [The same remark, of course, also 
applies to (32).] 

We are therefore always at liberty to add the same constant to all 
the p*(k)'s. Under some circumstances, this degree of freedom, together 
with the one provided by the introduction of k, can be used to reduce 
the dimension of B' and B : entire rows in B' and the corresponding 
columns in B can be set equal to zero without affecting the computation 
either of (31) or (32), or of the variance of the switch-count load. It 
is relevant to note here that this reduction would be largely illusory 
were it not for the fact that the normalized components of any of the 
characteristic vectors can be obtained without having to compute 
other components of that vector (see below). 

According to the result of Section 3.1, the covariance can always be 
cast in a form that involves only the first c components of the character- 
istic vectors. But when the input and departure rates for ^ n < c 
are linear in n, the covariance can also be expressed in terms of the 
last d — c + 1 components of these vectors. Indeed, the rates are 
then of the form 

X„ = \n + X', 

Hn = fin + m', n = 0, 1, • • •, c - 1, 

so that 

A„ - fin= (X - n)n + (X' - n'). 

Hence, with k = X — n and y = n' — \', (27) yields 

f0 if n - 0, 1, •••, c - 1, 

P«(X - m) = \ in - X)c - He + m' - ^' if n = c, 

[ (ji - \)c + h' - X' if n = c + 1, • • •, d. 

For the random (Poisson) and the quasi-random inputs, the p's take the 
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following simple form whenever the service time is exponential with 
mean 1. 

(i) Random input (X„ = a, >i = 0, 1, • • •) : 

f if n = 0, 1, • ••, e - 1, 

p*(-l) = \ -a if /( = c, 

v — a if n = c + 1, • • • , d. 

(ii) Quasi-random input [.V sources, X„ = (.V — n)\, 
n =0, 1, •••,-Y]: 

TO if n = 0, 1, •••,{ - 1, 

p n *[-(l + X)] = ! (c- -V)A if n = c, 

c + (c - N)\ if // = r + 1, • ••, (I. 

From the preceding developments, we see that the p's can be chosen 
in such a way that the number of components of the characteristic 
vectors needed to express R c is the smaller of the two integers c and 
d — C + 1. In particular, in the case of loss systems, only the (c +l)st 
component of each vector is needed. 

The parameters k and 7 can also be chosen so that only the first 
c + 1 components of the characteristic vectors actually enter in the 
expression of A',.. This will be (he case if we set k = n,-\ — A,- 1 and 
7 = c(mc-i - K-i). 

In Ref. 3, the derivation of the covariance function for loss systems 
[(/ = c, N e (t) = A r (0] with Poisson input and exponential service 
time makes use of the differential equations 

^-EiN(t)\N(0) = m] = - E[N(t)\N(0) = m] + a[l - p««(0], 

»; = 0, 1, • • •, c. 

These ecpiations appear here as that particular instance of (26) for 
which k = — 1, X„ = a, n — 0, 1, • •, c— 1, and/Z" = n, u = 1, • • •, c. 
Note also that now £„ f =o />„,„ (7) = 1 — p mc (t). But we stress that, in 
Ref. 3, the determination of the covariance relies on known recurrence 
relations between the so-called "sigma" functions (Ref. 13, pp. 129 and 
143 (T.); the more general problem considered in the present paper is 
not as readily amenable to such a treatment because of the greater 
complexity of the expressions that would now have to be used instead 
of the sigma functions. As we have seen, however, relatively simple 
formulas for R c can be obtained without extensive algebraic develop- 
ments as long as the underlying process is reversible. 
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IV. VARIANCE OF THE SWITCH-COUNT LOAD 

The variance of the switch-count load is now readily obtained. 
Depending on which expression we select for R c , we have either 

(i) Var L„(T) 

= nr* £ (ft- \k\)R c (k T ) 

k=-n 

= n- 2 -q' d -B' d -diag \o, ^^ (n - |fc|)e" |fc ", i= \, ■■■, dj 

•Brf-qJOO, (34) 

with 

q = (popl •■-, Pc-ipLi, 0, ■ • -, 0), or 

(ii) VarL„(r) = n-^KWJ-Bi 

■B^-qJW, (35) 
with k # r,-, j = 1, •••, d, or 

(Hi) VarLn(T) = ft- 2 {cr 2 c - [qSOOJ-B, 

•diag<» [0, (r,- - k)- 2 , i = 1, • • ■, d]-B rf .q rf | 

• £ (ft- |fc|)e«l*i r + n- 2 [q^(K)]'-B' d 

*=- n 

[n e r,|t|T "I 

0,Z_ n (n- \k\) {r _ _ K)i ,i= 1» ■■-,«*] 

•B,rq:W, (36) 

where k ^ r; for i ^ j". 
We now make use of the following identity (Ref. 3, p. 137) : 

* 1 — e~ 2nu 
£ (n - |fc|)e- 2|fc|u = n-cothw 5 csch 2 w. 

n = — k 

By means of this relation, (34) to (36) can also be written as 

[/ _ rr . \ i — e nTTi 
0, coth ^ -y J ^ — 

■csch* (^pYi= 1, ••,dJB d q d ; (34a) 
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VarL.CT) = ir 1 - [qj («)]'• B;, 

diag [O, ^-^ |coth (=p)- l -=$P 

k ^ r,-, ?'= 1, •••,<-/; (35a) 
and 

Var Z,«(T) 

= "- 1 {^-[q:W]'-B;-diag^ [0, {r t -K)-\ i=l, ■ ■ ■, dJ-B^qJW} 

•jcoth(^) - L^_p. csch 2^M + n -».[q:W7-Bi 

•cscW^p')! ,r=l, ••• > rfl-B tf -4i(«), ( 36a ) 

where k ?^ r,- for i 5^ ./. 

Let VarL 00 (7 1 ) = lim.^ M Var L n (T) be the variance of the load 
measurement obtained by continuous observation of the number of 
busy servers. If we replace r by T/n in (34a) to (36a) and then let n 
tend to infinity while keeping T fixed, we obtain the following formulas : 

Var L M (T) = - |<fc-Bi 

•diag[o,i(l + l ~Tr< T )> i= l > ■••»^]"B*-to (34b) 

VarL M (7 1 ) = - | [qJWJ-B; 

diag [°' r.-Cr,--.)' + ^T"" ) ' ' = lf ' ' '' f/ ] B " ,q2(<) ' 

k ?^ ?',, t — 1, • • •, d, (35b) 

and 

VarL„(T) = \\c\- [qSWl'-Bi 
•diag^[0, (r.- -j:)- 8 , »- 1, •■•■«*] 

•B^.qswj i(i + ^r^) - Ikgot-b; 

diag- [0, ^^ 2 (l + l -^f ) , * = 1, • • -, d] .B d .q; W , 

k ^ n for i ^ j. (36b) 
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We note that the formula for the variance of sums of dependent 
random variables makes it possible to compute the covariance between 
load measurements performed over distinct time intervals. Indeed, 
consider for instance a sequence of load measurements over the 
intervals (0, T], (T, 2T], (27\ 37], • • •. Let L™(T) be the switch- 
count load over the ith interval (i = 0, 1, • • •), S«(0 = » 2 ^«(0, and 
r« = Cov [L."»(T), I?(T)J Then 

Var £„<*+»[(* + 1)7] - (k + 1) Var5 n (T) 



t 



+ 2KH1 -i)n»r«(T) 



and 

r »(D = i^+Jl 2 Var L B(fc+l )[(* + 1)T] - Mjp Var L„(T) 

- ir,»(T). 

The preceding formulas may be used to determine the I\P*(T) re- 
currently. But the results of such computations shall be exact only if, 
for some choice of the time origin, all the scanning instants are multiples 
of r. 

We conclude this section with the remark that the variance formulas 
(34), (34a), and (34b) are valid for arbitrary reversible Markov pro- 
cesses with finite state spaces. 

V. NUMERICAL CONSIDERATIONS 

The exact variance formulas of the preceding section are very well 
suited to electronic computation and are easily programmed since, 
apart from straightforward evaluations of hyperbolic functions and 
simple products of matrices and vectors, they only involve the deter- 
mination of characteristic values and vectors for which powerful sub- 
routines are readily available. The fact that S d is symmetric and tri- 
diagonal (or reducible to tridiagonal form by an orthogonal similarity 
transformation) allows us to use the subprogram TQL2, which is par- 
ticularly efficient under the present circumstances (Ref. 11, pp. 227- 
240). Without going into details, we mention here only that this sub- 
program is based on the so-called QR-algorithm and relies on the con- 
struction of a sequence of symmetric tridiagonal matrices, S£ n) , 
n = 1, 2, • • •, unitarily similar to S d , which converges to diag [0, n, 
r>, ■•-, rj. At the nth iteration S£° is expressed as a product of an 
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orthogonal matrix Q d n) and a lower-triangular matrix V d n) : 

In Ref. 11, this decomposition is carried out by the Givens' triangular- 
ization in which the Q d n) 's are expressed as products of simple plane 
rotations (Ref. 7, pp. 239-240). 

The (?i -f- l)st iterate of S d is then given by 

whose unitary similarity to S d n) (and hence to S d ) follows from the 
relation 

UP-QP = {QPy-SP-QP. 

This method avoids the numerical difficulties frequently associated 
with the computation of the zeros of the characteristic polynomial. 
As shown in Ref. 11, p. 228, 



Qr = Rr 



R(D 
a » 



(37) 



where Ri" is of the form : 



1 






row i. 



(38) 



According to Ref. 11, p. 231, the matrix B' rf (whose columns are the 
characteristic vectors of S d ) is given "almost to working accuracy" by 



<?r 



QUI 
d j 



(39) 



where n (^30) is the number of iterations needed for the (numerical) 
symmetrization of S d . Taking (37) to (39) into account, it is then 
readily seen that the elements of B' d can be determined row by row 
which, as we have remarked earlier, is a desirable feature in the present 
context. 

Computations have been carried out for systems having as many as 
400 devices (and hence transition-rate matrices of order 401) to deter- 
mine the numerical accuracy of the approach described in the preceding 
sections. Checks were performed by comparing the value of Var Li(T) 
obtained by means of (35a) or (36a) with the corresponding <j 2 c , which 
can be calculated directly and independently from the equilibrium 
state probabilities. These two quantities, which are theoretically equal, 
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turned out in all cases to agree to at least 10 decimal places with the 
greatest difference occurring when d was largest. Hence, our procedure 
indeed yields very accurate results for the type of systems that are 
likely to occur in practice. But when d is large, the storage require- 
ments and the amount of computations become critical. It is therefore 
always important to select k and 7 in such a way as to minimize the 
number of B' rows that actually enter into the computations. (It 
follows from earlier remarks that this number, for proper choice of k 
and 7, never exceeds the integral part of (d + l)/2.) Further reduction 
can also be achieved by excluding the states whose probabilities of 
occurrences are so small that neglecting them will not materially affect 
the final results. In this connection, we make the following remarks. 
The variance of the switch-count load is perturbed by at most 

[p*00] 2 -P,-* 2 c 

if pj is set equal to zero in the particular formula used to evaluate 
Var L„(T). Hence, since 

Var L„(T) ^ a-Jn, 

we always have the following upper bound for the relative error, e J7 
induced by setting pj equal to zero : 

e> ^ [p*W] 2 -Pi-w, j = 0, 1, •••,(/. 

For a given relative accuracy of Var L„(T), these inequalities make it 
possible to determine ahead of time whether some components of the 
characteristic vectors can be "safely" eliminated from the computa- 
tions. In large systems, the gains achieved by such a reduction may be 
quite substantial, as either low occupation states [iV(0 small] and/or 
high occupation states [N(t) large] have then frequently very small 
probabilities of occurrences. 

Computations could be arranged to determine only those character- 
istic roots that are required to reach a given degree of accuracy 
[plus those needed to compute Bj*q2(«)D- This is rather readily 
achieved in loss systems with Poisson input and exponential service 
times since, in this case, the coefficients b] of 

in the variance formulas of Section 3.1 are then monotonically de- 
creasing as I Vi | increases : 

b\ < b 2 j if |r,| > \rj\, i,j = 1, ■ ■ ■ , d. 
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d = 80 
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ROOT INDEX 

Fig. 1 — Spectral measure of the carried load process. 

But, in general, one would encounter an additional difficulty, namely, 
that the b 2 t 's do not have the monotonicity property alluded to above 
and may actually fluctuate widely. This is illustrated in Fig. 1, where 
the roots are assumed to be indexed in order of increasing magnitude 
and the ordinates are the corresponding b^'s, normalized in such a way 
that max.- b\ = 1. 

The computations should be based on (35a)- (36a) or on (35b)- (36b) 
in the case of continuous measurements — as these formulas provide 
us with all the flexibility needed to cut down both storage space and 
computation time. When choosing between (35a) and (36a) or between 
(35b) and (3Gb), one should keep in mind that, for k close to r,-, the 
difference /', — k may not be determinable with enough precision to 
allow accurate computation of Var L U (T). This is shown in Table I 
where k = — 1 and n is the root of smallest positive absolute value. 
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Table I — Loss system, 80 servers, Poisson input, 
exponential service 



Offered 


1 + r, 


ffc* 


VarL,(T) 


Load in 
Erlangs 


Formula (36a) 


Formula (35a) 


10 
20 
30 
40 


1.09 X 10-' 3 
-9.90 X 10" 13 
-1.14 X 10~ 11 
-3.38 X 10-« 


10.000000 
20.000000 
30.000000 
39.999986 


10.000000 
20.000000 
30.000000 
39.999986 


0.034614 

0.013173 

20.891365 

39.999926 



(Note that the last two columns of this table should be equal and that 
errors of the same magnitude would arise if one were to use (18) of 
Ref. 3.) In all our computations, we have made use of (35a) and (36a) 
whenever \k — rt\ < 10 -4 for some i. This bound for \k — n\ is both 
large enough to ward off appreciable accuracy losses and small enough, 
under prevailing conditions, to be satisfied by only one root. 

VI. REMARKS ON INFINITE SYSTEMS 

It is known that infinite systems can be regarded as limits of finite 
ones, 12 and it is therefore of practical interest to have information 
concerning the spacing of the characteristic values as the dimension, 
d, of these approximating systems becomes large. Indeed, as d increases, 
computational difficulties may arise because of a lack of separation 
between these roots. Such problems would certainly come up sooner 
or later if the spectrum of A = lim^A,* happens to be dense over 
some interval as, for instance, in the case of a single-server queue with 
Poisson input, exponential service time, and unlimited waiting room 
(Ref. 12, pp. 365-366). Infinite systems with well-separated roots do, 
of course, also occur. As an example of this type, we mention the 
systems with an infinite number of servers, Poisson input, and ex- 
ponential service which often provide useful idealizations. (In this 
case, as is well known, the nonvanishing characteristic roots are the 
negative integers, -1, -2, -3, • • •.) Other examples of systems with 
discrete spectra are given in Ref. 12, where sufficient conditions for 
this to occur are discussed in some details; but in all these instances 
the X„ and n» both increase as n" for some v > 0. This condition is un- 
likely to be satisfied in queuing systems ; generally, in this particular 
area, the arrival and the departure rates remain bounded : 



^ X„ ^ A < oo, ^ /m £ M < », n =• 0, 1, 



(40) 



As briefly described below, these inequalities imply the existence of 
definite bounds for the spectrum of A. 
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Consider an infinite system, and let A be its (infinite) transition-rate 
matrix. Let A d be the matrix obtained by retaining only the elements 
belonging to the first (d + 1) rows and columns of A and then setting 
\a equal to 0. Let rdo( = 0) > rai > • • • be the characteristic roots of 
Ad. Then, under conditions (40) it can be shown that, for any k ^ : 

(i) | rrf* | < A + M for d sufficiently large, 

(it) kd. d - t | < 2 (A + M) for d ^ k. 

Either of these two inequalities implies that the characteristic roots 
do not remain separated as d —* °o whenever (40) is satisfied. Under 
the more stringent requirements that (40) holds and that 

lim X„ = A, lini/U" = M, 

n-»oo n-*oo 

more precise statements can be made, namely, that, for all k's and d's, 

\r dk \ < (VA + VA?) 2 
and that the spectrum of A always comprises a closed interval, viz., 

n = [_ (VX + Vm) 2 , - (Va - Va?) 2 ]. 

(In addition to £2, the spectrum of A may also include a finite number of 
roots in [— (VX — Vl?) 2 , 0].) But it turns out (as will be shown 
elsewhere) that, as d increases, the characteristic roots of Ad fill fi 
rather "evenly"; furthermore, for practical accuracy levels, large 
values of d are needed only when the length of Q, tends to be relatively 
large (a circumstance corroborated by extensive computations). 
Hence, within the present framework, it appears that root-spacing is 
not likely to be critical except in the improbable event that extreme 
precision is required. 
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